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Over the past few years new physics methods and algorithms as well as the latest snpercomputers 
have enabled the stndy of the QCD thermodynamic phase transition nsing lattice gange theory nn- 
merical simulations with unprecedented control over systematic errors. This is largely a consequence 
of the ability to perform continuum extrapolations with physical quark masses. Here we review re¬ 
cent progress in lattice QCD thermodynamics, focussing mainly on results that benefit from the 
use of physical quark masses: the crossover temperature, the equation of state, and fluctuations of 
the quark number susceptibilities. In addition, we place a special emphasis on calculations that are 
directly relevant to the study of relativistic heavy ion collisions at RHIC and the LHC. 
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I. INTRODUCTION 

Quantum Chromodynamics (QCD), the prevailing theory of the strong interaction in nuclear physics, is one of the 
most successful and most challenging theories in physics. At low temperatures, chiral symmetry breaking gives rise 
to most of the mass in the visible universe, and at high temperatures it predicts the existence of a plasma of quarks 
and gluons BE]: a condition achieved in the first microsecond after the big bang as well as in nuclear collisions at 
the Relativistic Heavy Ion Collider (RHIC) and the Large Hadron Collider (LHC). However, calculations at finite 
temperatures in the vicinity of the transition between hadronic matter and quark gluon plasma require computationally 
challenging techniques in lattice gauge theory Bill. Steady progress in this field relies upon continuing advances in 
both algorithms and computing. Recent improvements in fermion actions and computing capabilities have enabled 
lattice gauge calculations to perform the first reliable continuum extrapolations with physical quark masses. The 
QCD transition observed in heavy ion collisions is now firmly established as a crossover at zero baryon density, 
and several groups have produced continuum extrapolations for the crossover temperature and equation of state 
(EoS) with physical quark masses. These calculations are consistent within their respective uncertainties, which now 
have well defined statistical and systematic contributions. Furthermore, access to high statistics data from heavy 
ion collisions coupled with new analysis techniques and improvements in hydrodynamic modeling tools have greatly 
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enhanced the connection between the lattice EoS and experimental data. Similar achievements are expected soon 
regarding calculations of heavy quark color transport and screening. Finally, new insights into the role of conserved 
charge fluctuations on the lattice and in heavy ion collisions are providing additional avenues for direct comparisons 
between calculations and data. 

In this review, we summarize recent results in calculating the basic thermodynamic properties of high temperature 
QCD, with special emphasis given to calculations with physical quark masses. These include calculations of the 
crossover temperature, the equation of state, and quark number susceptibilities. We show how these results play a 
crucial role in understanding the properties of the quark gluon plasma created in heavy ion collisions. We also provide 
a brief review of the current status of heavy quark color screening and transport, for which physical quark mass 
calculations are not yet in reach, and we discuss implications for future progress. In Section |ll] we describe the gluon 
and fermion actions used in recent QCD thermodynamics calculations and recent results are presented in Section [III| 
In Section |IV| we summarize recent results for color screening and transport, and in Section [V] we offer conclusions 
and discuss future prospects. 

This review is aimed primarily at students entering the field and scientists outside the field who are interested in 
recent results and future possibilities for lattice QCD calculations in thermodynamics. With this goal in mind, we 
anticipate that experimentalists and non-lattice gauge theorists in relativistic heavy ion physics will find this review 
especially useful. 


II. LATTICE QCD ACTIONS AND PROPERTIES 

The study of QCD thermodynamics begins with the partition function, expressed as a path integral over the classical 
Euclidean action separated into fermionic and gluonic components, 

Z{T,V)= J DA^D^D^eM-S^co) ( 1 ) 
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for temperature, T, with the sum over ny different quark flavors and implicit sums over double-indices for color 
degrees of freedom of gluons (a = 1, ..N^ — 1) and quarks (a, (3 = 1, ..W)- In the discretized (lattice) formulation, the 
integrals are replaced by sums over fermion fields occupying each lattice site and compactified gauge variables on the 
links connecting fields on neighboring sites. The sums are performed over spatial and W temporal steps, related 
to the temperature, T = l/iaNr), where a is the lattice spacing. In many thermodynamic calculations it turns out 
that an aspect ratio for N^j/Nt of 3 or 4 is already close to the infinite volume, thermodynamic limit. Therefore, a 
given lattice calculation is referred to by specifying the number of temporal steps, W. The relative contributions of 
fermionic and gluonic terms vary with the observable and temperature range. For the trace anomaly from which the 
EoS is derived, the gluonic term contributes 80% near the crossover, increasing to 90% by 400 MeV. However, the 
physical contributions from quarks and gluons do not separate so cleanly. For example, the fermion term vanishes for 
massless quarks although they certainly contribute to the EoS. 




FIG. 1: Left: Diagrammatic representation of terms in the improved lattice gauge actions. The three vector diagrams on the 
left depict gluonic links connecting neighboring space-time lattice sites. The left most represents the unimproved Wilson action, 
referred to as a “plaquette”, and those to the right represent increasing levels of improvement. The first diagram on the right 
(Cl) represents the unimproved staggered fermion action and the remaining diagrams (Ca-Cjv) are included as improvements. 
See text for further details. 


Over time, additional terms have been added to the action to remove errors of various orders in the lattice spacing, 
thereby improving convergence to the continuum. Examples of gluon terms are shown in the left three diagrams of 
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Figurej^ Closed paths (Wilson loops) generate different approximations to the continuum gauge action term . 

The plaquette on the left (/3) generates the unimproved Wilson action. Symanzik-improved gluon actions include 
additional terms in the action that remove errors of 0{a^) at tree level [S]. This term is represented by the rectangle 
shape {Prt)- With the addition of the “chair” {j3ch), diagram to the right, sometimes called the “parallelogram”, it is 
possible to eliminate, errors of 0{asa?)^ resulting in the one-loop tadpole-improved action of Liischer and Weisz [S]. 
Nearly all contemporary thermodynamic calculations employ some version of Symanzik improvement for the gluon 
action. 

Most of the recent algorithm improvements in lattice gauge theory have been devoted to improvements to the fermion 
action. A naive discretization of the Dirac equation for fermions introduces additional minima in the dispersion 
relation at momentum component tt/o, effectively doubling the number of fermions for every dimension, i.e. 2^ 
for each flavor. There are several approaches to the fermion doubling problem. The action introduced by Wilson [5] 
includes an additional term that imparts a very large mass to all but one of the fermion eigenstates, but at the expense 
of explicitly breaking chiral symmetry. Because of the important role of chiral symmetry in finite temperature QCD, 
Wilson fermions are used less often in thermodynamics, where actions that preserve at least some aspects of chiral 
symmetry are preferred. 


A. Staggered Fermions: p4, asqtad, stout, HISQ 



FIG. 2: The root-mean-square (RMS) mass of members of the pion taste multiplet for asqtad, stout, and HISQ/tree actions 
as a function of lattice spacing. These actions employ different improvements to the staggered fermion action. The HISQ/tree 
action has the lowest RMS mass over this range. 


The staggered fermion formulation exploits an exact lattice symmetry to reduce the number of unwanted doublers in 
the naive action from sixteen to four per flavor. In effect, the Dirac spinor components are interleaved on alternating 
lattice sites. For each flavor this introduces an additional, unphysical quantum number referred to as “taste”. A 
part of the chiral symmetry on each site is preserved, and full chiral symmetry is restored in the continuum limit. 
In this limit, the taste degrees of freedom are 5'C/(4)-symmetric, and the correct number of degrees of freedom is 
restored by taking the fourth root of the fermion determinant for each flavor. The difficulty arises for finite lattice 
spacing, where the taste symmetry is broken. The impact of taste splitting appears in the light meson spectrum. 
Because a remnant of chiral symmetry is preserved, the lightest pion behaves as the traditional Goldstone boson 
whose mass vanishes in the limit of zero quark mass. The remaining 15 members of each taste multiplet have seven 
partially degenerate masses, but in the limit of zero lattice spacing, all taste-symmetry-violating splittings vanish. 
This situation is illustrated in Figurej^which shows the root-mean-square (RMS) value of pseudo-scalar meson masses 
(would be pions) for several different staggered fermion actions to be explained below. For physical values of the light 
quark masses, the different masses from taste splitting converge towards the pion mass in the continuum limit. 
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Improved fermion formulations include additional terms to reduce lattice artifacts (“lattice cutoff effects”), one 
manifestation of which are the taste splittings. The Symanzik improvement strategy adds terms to the fermion 
action so that, in the continuum limit, it still reduces to the correct Dirac action, but errors of 0{a^) or higher are 
eliminated. This is the philosophy of the asqtad [7] and p4 [5] actions, which have been used extensively for QCD 
thermodynamics. Diagrammatic examples of the asqtad improved fermion action are labeled {Ci-Cn) in Figure 
The terms depict fermion hopping between next-neighbor and third-neighbor sites. The paths of gluonic links define 
the gauge connection between fermion fields on neighboring sites. The single Ci link on the left yields the unimproved 
staggered action with discretization errors of 0{a^). The five other path shapes are needed to eliminate those errors 
completely, leaving errors of 0{a‘^) and ©(OsO^). The p4 action is similar, but replaces the three-link term with a 
“knights move” term and includes other small changes [8]. 

Further improvements are possible. The asqtad and p4 modifications reduce taste symmetry breaking by suppress¬ 
ing hard gluon exchanges that cause transitions between tastes. They do this by smoothing the local gluon field 
experienced by the fermion. The highly improved staggered quark (HISQ) action has a still higher level of smoothing. 
In a different approach, rather than attempting to achieve an exact cancellation of terms of 0{a^), the Budapest- 
Wuppertal group smoothed the gluon field “stout smoothing” [3] before coupling it to otherwise unimproved staggered 
fermions. 

Another important measure of improvement, particularly relevant at higher energies where quark degrees of freedom 
predominate, is the extent to which the quark dispersion relation is accurately represented. The asqtad, p4, and HISQ 
actions include terms that eliminate O(a^) corrections to the quark dispersion relation. Thus, they are expected to 
perform well at high energies. The stout action does not include these additional terms, opting instead to perform 
calculations with less improvement at smaller lattice spacings. However, in the continuum limit all cutoff effects 
disappear and the various staggered actions should all agree. 


B. Chiral Fermions: domain wall and overlap 



FIG. 3: Illustration of fifth dimension of the domain wall fermion action, in which the left (4 'l) and right-handed (4'i?) chiral 
states are exponentially bound to opposing walls separated by a distance of Ls lattice sites. The red region indicates the small 
overalp that occurs at finite Ls and is responsible for mixing the chiral states and producing a residual mass rrires. An explicit 
input mass m/ that directly couples the two chiral states is introduced to provide explicit control over the mass. The “height” 
of the domain wall is denoted by mo. 


To enable a fully chiral treatment of fermions on the lattice one must be able to independently rotate the two 
chiral states even at finite lattice spacing. Domain wall fermions (DWF) achieve this by binding the two chiral states 
to opposing domain walls in a fictitious fifth dimension [T(1IH2j . Figure illustrates how the left and right handed 
chiral states are exponentially localized on opposing walls, thereby achieving chiral symmetry even at finite lattice 
spacing. Chiral symmetry is now restored separately from Lorentz symmetry. The former is fully recovered at the 
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Ls —>■ oo limit {Lg is the size of the fifth dimension) while the other at the a —>■ 0 limit. This eliminates the taste 
splitting that exists for the staggered fermion formulations and instead gives rise to exactly 3 light pions, enabling the 
study of subtle features of chiral symmetry restoration. However, the extra dimension brings with it an additional 
computational cost which increases linearly with Lg. 

As shown in Figurej^the left and right chiral wave functions have some overlap and therefore break chiral symmetry. 
Because they are expected to be exponentially localized on opposite walls, this overlap is expected to be small. 
However, this overlap induces an additive mass, called residual mass mres, to the input quark mass ruf, and as a 
result TOeff = mf+rrijies- The primary challenge in applying the DWF action to lattice QCD is in achieving small values 
of rripes while controlling the overall computational cost that grows with Lg. This is especially true for thermodynamic 
calculations, which require zero-temperature calculations on large volumes in order to eliminate divergent vacuum 
contributions. In the domain wall formulation, the gauge links are replicated on each four-dimensional slice while the 
links in the fictitious fifth dimension are all set to the unit matrix. This allows the construction of a transfer matrix 
T and corresponding Hamiltonian H 4 . Such a Hamiltonian framework is called the overlap formalism m- Variants 
of the overlap formalism have been developed that make it suitable for numerical simulations [141116) and have been 
used in various Lattice QCD studies. These variations have different technical properties and problems than DWF 
but are exactly identical at the infinite Lg limit. 

In the DWF approach, the localization of the two chiral components is controlled by the eigenvalue spectrum 
of the transfer matrix Hamiltonian (see for example [17]). Since the states on the opposite walls are related by 
exp(— Lg X iJ 4 (mo)) their localization depends on the eigenvalue spectrum of Hi^mo). The lattice QCD simulation 
generates an ensemble of gauge field configurations, and the eigenvalue spectrum of 774(mg) will be different for each 
configuration. As a result the localization and m-res will vary from configuration to configuration. For any gauge field 
configuration, H^^mo) has the same number of positive (?t,_|_) and negative (n_) eigenvalues for toq < 0. However, 
as the height of the domain wall, mg, is increased above zero some eigenvalue of 774(mg) may cross zero and change 
sign. Then the quantity (n+ — n_), which is the index of the Dirac operator, would not be zero just after the crossing 
occurs. It has been shown that the number and direction of crossings is directly related to the number of instantons 
and anti-instantons present in the gauge configuration and that (n_|_ — n_) is equal, in a statistical sense, to 
the net (global) topological charge of the gauge field configuration [TSj. This is of particular importance to faithfull 
lattice studies of the 77(1)Axial symmetry breaking in QCD. However, it follows that a configuration that is close to a 
topology change will have a near zero 774(mg) eigenvalue and therefore poor localization and large mres- 

This limitation has recently been addressed by applying a Boltzmann weight |77|(mg) + e'^\ to the action. The 
771 (mg) term provides a better localization and smaller mres |19j . and the e term allows one to control the topological 
charge crossings. This method was first used for the study of QCD thermodynamics in |2()j and has been named 
DSDR (dislocation suppression determinant ratio) [21]. Because of the created gap one is able to reach small enough 
mres to produce three degenerate pions at their physical value of about 140 MeV and faithfully simulate the QCD 
thermal transition |52|. The DWF algorithm continues to benefit from algorithm improvements; see [23] for a recent 
review. 


C. Temperature Scales and Cutoff Effects 

The lattice calculations are performed in dimensionless units, and converting the coupling parameter to physically 
meaningful temperature scale requires matching either directly or indirectly to an experimentally measured quantity. 
Two common examples are the rg scale set by the T (2S-1S) mass splitting, and the meson decay constants, 

/k- The former is determined by setting the derivative of the static quark potential = 1-65 for r = rg. A 

second value ri for which this quantity is equal to unity is also used. This scale is more suitable for calculations with 
smaller lattice spacings as cut-off effects start becoming small at this distance and ri is statistically more accurate 
than rg. The decay constants are more influenced by cutoff effects that also affect the trace anomaly in the vicinity 
of the transition temperature where this quantity is dominated by contribution from hadron resonances. Therefore, 
when using the fx scale, thermodynamic quantities are more similar for different lattice spacings and the continuum 
extrapolations are less severe. Recently a new Wilson flow scale, wq has been proposed, which makes use of the D 
mass for setting the scale |24|- At this time, with lattice calculations performed much closer to the continuum limit, 
the choice of scale parameter is less controversial than it once was. Systematic errors associated with the scale setting 
are now only 1-2%, and many results are cross-checked using more than one scale setting. 

The term “cutoff effect” is used to describe any error associated with a particular lattice discretization scheme 
used in a calculation. We have discussed the taste symmetry breaking of the staggered fermion action, but there 
are additional corrections depending upon the observable being calculated. Although all such artifacts are expected 
to disappear in the continuum limit, they can lead to significant deviations from the scaling approximations used 
to perform the continuum extrapolations. Knowing when the scaling region has been reached is perhaps the central 
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challenge in any lattice calculation, and it is only truly evident in hindsight, when an additional calculation at smaller 
lattice spacing produces no significant change in the extrapolation. 


III. LATTICE QCD THERMODYNAMICS RESULTS WITH PHYSICAL QUARK MASSES 

A. Transition Order and Temperature 



FIG. 4: Order of the QCD phase transition vs. strange and light quark mass values. 


The location of the QCD “transition temperature” has received much attention in the heavy ion community, but a 
single temperature can be defined precisely only for a true phase transition, in which the singularity in the partition 
function extends to all observables. However, understanding the nature of the transition and defining its corresponding 
temperature in a meaningful way have important implications for heavy ion collisions and our overall understanding 
of QCD thermodynamics. 

Figure illustrates the nature of the phase transition as a function of quark mass [551 US], in which the phase 
boundaries are motivated by lattice calculations and symmetry arguments. It illustrates the importance of having 
lattice calculations with physical quark masses, which now confirm that the transition is a crossover at the physical 
point of the diagram. 

The most direct method to probe the nature of the QCD phase transition is to study the derivatives of the log of 
the partition function with respect to the light quark mass: the chiral condensate, which becomes the order 

parameter in the chiral (massless) limit, and the chiral susceptibility, Xtjjiin 


TdlnZ Td^lnZ 


(3) 


For a true phase transition the chiral susceptibility becomes narrower and the peak height increases with increasing 
lattice volume. The peak height increases linearly with volume for a first order transition and grows with critical 
exponents for a second order transition. Figure (left) shows calculations with the stout-link improved fermion 
action m- These calculations with physical quark masses show no volume dependence, but calculations with staggered 
fermions leave open the question of a potential systematic error associated with the lack of a true chiral symmetry on 
the lattice. 

This question has recently been answered with a calculation using domain wall fermions, shown in the right panel 
of Figure Lattices with the same lattice spacing and spatial dimensions of 32^ and 64^ have identical values of the 
disconnected chiral susceptibility. 






LQCD Thermodynamics with Physical Quark Masses 


7 


A/t = 6 




T (MeV) 


FIG. 5: Chiral susceptibility for several volumes for the stout staggered action for Nt = 6 with physical quark mass (left), and 
the domain wall fermion action for Nt = 8 (right) for rriT values of 140 and 200 MeV. Comparisons to the HISQ action with 
Nt = 12 are also shown. There is no evident change in peak height with increasing volume for the same pion mass. 


As noted, the presence of a crossover transition for physical quark mass values complicates the definition of a 
transition temperature, but many of the thermodynamic observables that develop singularities in the chiral limit may 
retain some remnant of the transition in a steep drop or inflection point in the crossover region, corresponding to the 
peak in the chiral susceptibility seen in Figurej^ These characteristics are used to define a pseudo-critical temperature 
(Tpc). 



T[MeV] T [MeV] 

FIG. 6: The renormalized chiral condensate (left) and the subtracted renormalized chiral condensate (right) for the stout action 
for Nt = 8,10,12,16 and the continuum extrapolation with physical quark masses. 


When working with the chiral condensate, it is common to remove lattice artifacts through subtraction and nor- 
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malization. The renormalized chiral condensate, {ilitp)^ and the subtracted 

= [W)i.r - W)i,o] 

W)i.T - 

W)/,0 - ’ 

where the I and s subscripts refer to the light and strange quark condensates and the T and 0 subscripts refer to 
the finite and zero temperature values, respectively. Calculations with the stout action for the renormalized chiral 
condensate are shown in Figure!^ with the renormalized chiral condensate shown on the left, and the subtracted chiral 
condensate shown on the right [28]. Results are shown for calculations with physical quark masses for Nj- = 8,10,12,16 
and the continuum using the Jk scale. Fits to the inflection point in {'4’4’) R yield Tp^ = 155 ± 2 ± 3 MeV, and similar 
fits to 6 i^s yield a value of Tpc = 157±3±3 MeV. One can also obtain Tpc by fits to the peak in the chiral susceptibility. 
In this case the authors fit to obtaining a value of Tpc = 147 ± 2 ± 3. 

The HotQCD collaboration has produced a similar result, Tpc = 154 ± 8 ± 1 MeV, using a method that relates 
the pseudo-critical behavior of the chiral susceptibility more directly to universal properties of the critical behavior 
of the true phase transition in the chiral limit. This is accomplished through the use of 0{N) scaling relations, which 
enable one to parameterize the behavior of the chiral condensate in the vicinity of the phase transition. The scaling 
relations are referred to as 0{N) because observations of both 0(4) and 0(2) scaling are possible. QCD belongs to 
the 0(4) universality class in the limit of vanishing light quark masses and a sufficiently large strange quark mass, 
but the staggered fermion action retains only an 0(2) global symmetry. The universal scaling relations and critical 
exponents for three dimensional, 0{N) symmetric models are well known and were first exploited for a discussion of 
the QCD phase diagram by Pisarski and Wilczek [29] . The parametrization of the scaling functions in a form suitable 
for QCD applications (SO] has been refined in recent years |3T| and has been previously applied to the p4 action [5^ . 

Because this approach is relatively new, and establishes a more direct link between Tpc and Tc in the chiral limit, 
we summarize the essential features of the scaling analysis used to perform a simultaneous fit to the asqtad and 
HISQ/tree continuum extrapolations for mi/rrig = 1/27. More details can be found in |33j and references therein. 

To isolate the scaling component, one separates the singular part of the 0{N) phase transition from the regular 
part. 


chiral condensate, Si^s, are defined as, 

(4) 

(5) 
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T 

--InZ. 



1 + 1/5 

y^singular ( T /regular (7^5 


( 6 ) 


where the singular part of the free energy density is expressed in terms of the scaling variable, z, which is defined 
in terms of dimensionless couplings related to the temperature, 7 = ^ , and the ratio of light to strange quark 

masses, h = ^ ^. The h parameter contains the quark mass dependence and plays the role of a symmetry breaking 
magnetic field. 


2 = 


I T-Tc 

zo Tc \Rni) 


(7) 


where /3 and 5 are the critical exponents. The chiral condensate can then be expressed as a function of the temperature 
and quark mass. 
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in which the regular part is expressed as a Taylor expansion to first order in t, and the coefficients oq and oi are 
determined from a fit. The parameterizations of the scaling functions for 0(2) and 0(4) have been obtained from 
calculations with three dimensional scalar 0{N) models |3T|. The results from fits to the HISQ/tree chiral condensate 
for Nr = 8 are shown in the left panel of Figure]^ and the corresponding functional forms for the chiral susceptibility 
are shown on the right panel. 

The quark mass dependence of Tpc is then given by. 


Tpc (h) = Tc 


1 + ^ 
Zo 


Qi^o 

y 2ApZpZo^-l^ 




-i/a-nZ/sA 


(9) 


where Zp and Ap are the peak position and curvature of the scaling function fxiz) = fx{zp) + Ap{z — Zp)'^, which can 
be approximated by a quadratic polynomial in the vicinity of the peak. This scaling analysis was performed for two 
lattice spacings for the asqtad action and three lattice spacings for HISQ/tree. These values were then incorporated 
into a simultaneous quadratic fit in V/, shown in Figure]^ 
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FIG. 7: 0(4) scaling fits and data for the chiral condensate Mb for the HISQ/tree action for Nt = 8 (left) and the corresponding 
values of the chiral susceptibility Xm,i (right). The difference between the peaks of the solid red and dotted black curves indicate 
how the regular part of the transition temperature changes with quark mass. 



FIG. 8: Simultaneous quartic continuum extrapolations of the transition temperature for 0(4) scaling for both asqtad and hisq 
actions [^. From this analysis, the pseudo-critical temperature for physical light and strange quark masses was determined 
to be 154 ± 8 ± 1 MeV. 


B. The f/(l)Axiai Symmetry 

The massless QCD Lagrangian also possesses an axial, 17^(1), symmetry. This symmetry is broken due to quantum 
fluctuations, giving rise to non-conservation of axial current [3T1 and leading to the explicit breaking of global 
C/a( 1) symmetry through the presence of topologically non-trivial gauge field configurations, namely the instantons 
[36]. With increasing temperature the density of instantons reduces drastically due to the color Debye screening [37] 
and eventually the 17 a( 1) symmetry becomes exact in the T ^ oo limit. Although, the 17 a( 1) symmetry is not an exact 
symmetry of QCD, the magnitude of its breaking around the chiral crossover temperature is expected to influence the 
nature of the chiral phase transition. As mentioned before, for two massless flavors and in presence of significant Ua{^) 
breaking the chiral transition in QCD is expected to be in the 3-dimensional 0(4) universality class |29l|38|. However, 
if the 1/^(1) symmetry breaking becomes negligible near the chiral transition temperature then the chiral transition in 
QCD can be of first order or of second order with the symmetry breaking pattern 0^(2) x Ur{2) —>• S'17y(2) [5^110] . 
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Therefore, determining the fate of [/a( 1) breaking at high temperatures provides a more complete understanding of 
the chiral transition. 

The issue of the axial anomaly for rooted staggered fermions is subtle and the correct anomaly may emerge only 
in the continuum limit SDIU]. On the other hand, the anomalous symmetry within the DWF formulation is more 
straightforward [43j : Ua{^) is broken by the same topologically non-trivial configurations as in the continuum and 
explicit lattice artifacts appear only at order Thus, the DWF action is a natural candidate to investigate the 

temperature dependence of D^(l) breaking in QCD. 

For two massless flavors the pion and the isovector scalar 5 (oi) mesons transform into each other via an Ua{^) 
rotation and the presence of an exact Ua{^) will render these mesons states degenerate. Thus, the difference of the 
integrated two-point correlation functions of pion and 5 meson, 

Xtt - Xs = j d‘^x[{-K+{x)TT~{Qi)) - {5+{x)5~{Q))] , (10) 


can be used as a measure of the Ua{^) breaking [33]. If the CfA(l) symmetry is exact then this quantity will vanish, 
and for small light quark masses, mp the non-vanishing corrections will be of the order of mf. This particular measure 
of 17 a( 1) breaking has been extensively studied using the DWF action as a function of temperature for different light 
quark masses and for several volumes [4^47] . As shown in Figure]^ (left), these DWF calculations clearly show that 
Xtt ~ Xs does not vanish around the chiral crossover temperature Tp^ = 154(9) MeV and remains independent of quark 
mass for T > 165 MeV, indicating that Ua{^) may remain broken at these temperatures even in the chiral limit. 
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FIG. 9: (Left) The figure shows that the Ua{X) breaking measure x-^ ~ Xs remains non-vanishing for T > Tpc = 154(9) MeV 
and becomes independent of quark mass for T > 165 MeV |47| . (Right) Infrared eigenvalue density of Dirac fermions that give 
rise to the Ua{1) breaking observed in x-n- ~ Xs |46| . 


The 1/^(1) breaking is intimately related to the topology of the gauge fields and consequently to the infrared modes 
of the Dirac fermions. In the limit of infinite volume, both the chiral order parameter, (i/'V'); and the Ua(X) breaking 
measure, ~ Xs^ can be written in terms of the eigenvalue density, p(A), of the Dirac fermions 



d\ 


2mip{X) 

+m’f 


and 


Xtt-XS = 



dX 


4mfp{X) 

(X^ + mf)^ 


( 11 ) 


In the limit m; —>■ 0 and for T > Tpc, ('(/''(/') vanishes but Xtt ~ Xs remains non-zero. This leads to an intriguing 
question: what is the form of p{X) that satisfy both these requirements? An answer to this question naturally leads to 
the underlying non-perturbative mechanism of axial symmetry breaking. As shown in Figure]^ (right), DWF studies 
suggest (35] 136] that an eigenvalue density of the form p{X) ~ mfS{X) can largely account for the C7 a( 1) breaking 
observed in Xir — Xs- Such an eigenvalue density naturally arises within a dilute instanton gas— a gas of widely 
separated, weakly interacting small instantons and anti-instantons. A more recent study using overlap fermions, 
possessing even better chiral properties and an exact index theorem, have put this underlying mechanism of C7 a( 1) 
breaking at high temperatures in a much firmer footing |48] . 
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C. QCD Equation of State 


The last few years have produced two full continuum extrapolated equation of state calculations at or near the 
physical pion mass: one each from the Wuppertal-Budapest [49j and HotQCD [50] collaborations. The results, 
shown in Figure]^ are consistent within systematic errors in the region from 130-370 MeV. These calculations were 
performed with the stout and HISQ/tree actions described in Section |n[ The overall consistency between the two 
results is remarkable given the differences in the staggered fermion actions and analysis methods used to extract the 
EoS in the continuum limit. 

The equation of state is derived exclusively from the trace of the energy momentum tensor, 0^^, referred as the 
trace anomaly or interaction measure because it serves to measure deviations from the conformal equation of state, 
in which the energy density e is equivalent to three times the pressure P, 

e- 3P r d ( P 



On the lattice, the trace anomaly is defined by the derivative of the log of the partition function with respect to the 
lattice spacing, ~ f dinf ’ evaluated separately for the gluonic and fermionic operators, 


e-3P_0^'"(T) ^ 0>^/{T) 

y4 2^4 y4 ’ 

^!^=Rp[{sG),-{sG)T\Nt, 


0 ^/{T) 


Pms ((V’V’)s.o - ('0V')s.T)]^r. 


(13) 

(14) 


(15) 


where (sg) is the expectation of the gauge action, Rp = —a^ is the nonperturbative beta function, and Rm is the 
mass renormalization function. The gluonic component is calculated from the difference in the expectation values 
for the action at zero and finite temperatures, whereas the fermionic component is derived from the difference of 
chiral condensates evaluated at zero and finite temperature in the light and strange quark sectors, respectively. Thus 
for the EoS calculation, each finite temperature calculation requires an analogous computationally expensive zero- 
temperature calculation on a large lattice (i.e. 64^) with equivalent coupling. This is one reason why EoS calculations 
require large, sustained computing allocations and often take more than one year to complete. As noted previously, 
the division between gluonic and fermionic components is not strict, but the gluon contribution tends to dominate 
both the signal and errors, especially at higher temperatures. Further details on the EoS the determination of Rp 
and Rm are given in 0150]. 


From Eq. 12 it follows that the pressure is determined by integration over the trace anomaly. 


P{T) ^ ^ 

pi pi 


'To 


J'/5 


(16) 


and is therefore sensitive to the value used for the pressure at the lowest temperature. Calculations at temperatures 
in the low temperature region (significantly below Tpc) are prohibitively expensive and in this region the hadron 
resonance gas is expected to provide accurate estimates. Furthermore, a smooth transition to the hadron resonance 
gas EoS is required in order to match hydrodynamic outputs to hadronic cascade codes when modeling heavy ion 
collisions. The HotQCD collaboration uses the HRG EoS at 130 MeV as a matching condition for the continuum 
extrapolation, whereas the Wuppertal-Budapest collaboration performs an integration over quark mass to set the 
normalization of the pressure [51]. The pressure values calculated with each action/method differ by less than 10% 
over the full range of temperatures, well within their combined errors, and both actions achieve a smooth transition 
to the HRG EoS just below the transition. 

Both calculations were performed along the lines of constant physics (TCP), in which the quark masses are held 
constant at physical values. The Wuppertal-Budapest collaboration used the ratio of the pion decay constant to pion 
mass, Jk/M-k for this, whereas the HotQCD used the mass of the fictitious r]gs meson, . The 

overall temperature scale is set by Wuppertal-Budapest using /k, whereas the static quark potential and first derivative 
Top are used by HotQCD. However, both collaborations also calculate the scale with wq and obtain consistent results 
that are also incorporated into their respective systematic error estimates. 
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The two results also differ in the treatment of the continuum extrapolation and the estimate of systematic errors. 
The Wuppertal-Budapest continuum extrapolations are performed on lattices with Nt=6 , 8, and 10 with Nr=12 
included for three values of the temperature. Above 350 MeV, only the Nt=6 and 8 were used. The continuum 
extrapolation is performed on spline fits to the data, and the extrapolation is quadratic in the lattice spacing. 

Continuum extrapolations for the HISQ/tree action were performed on lattices with Nt=8 , 10, and 12 using 
simultaneous quadratic fit to splines in which the spline knot locations were included in the overall minimization 
procedure. The uncertainties were estimated by fitting 20k samples in which the lattice calculations were allowed 
to vary within normal errors. The match to the hadron resonance gas was achieved by sampling the HRG value at 
130 MeV allowing for 10% variation in the value and fixed slope at that point. Both collaborations have produced 
parameterizations of their EoS calculations for insertion into hydrodynamics models of heavy ion collisions. 



FIG. 10: Continuum extrapolation of the interaction measure, energy density and pressure for the HISQ/tree and stout 
actions (left), and the corresponding speed of sound for HISQ/tree and stout compared to the Hadron Resonance Gas at low 
temperatures. 


Final results for the stout and HISQ/tree action trace anomaly, energy density, and pressure are shown in the left 
panel of Figure The right panel shows the square of the speed of sound and a comparison to the HRG calculation 
at low temperature. The equation of state is an essential input for accurate modeling of heavy ion collisions with 
hydrodynamic simulations. The equation of state is used to convert the initial Glauber or glasma density profile to an 
initial temperature or entropy profile |52j . Thereafter only the speed of sound enters into the hydrodynamic calculation. 
Before lattice calculations were able to provide smooth parameterizations to the modeling community, an over-reliance 
on simple formulas, such as the bag model equation of state with first order phase transition added to the challenge 
of correctly modeling the space-time distributions at freeze-out, as measured by femtoscopic correlations |53l I54j . It 
was only through the simultaneous adoption of a lattice-inspired EoS and second order viscous terms that provided 
the modeling community with the tools needed to successfully model the evolution of a heavy ion collision [55] . 

The question at this time is whether current uncertainties are sufficient for current and future modeling needs, or 
whether additional refinements are needed (at significant computational expense). The answer to this question is not 
yet rigorously known, and depends upon current multiparameter sensitivity studies that are just beginning [5M58] . 
The prevailing consensus is that current uncertainties in the lattice EoS are sufficient and that further refinements will 
not greatly elucidate the dominant uncertainty in the understanding and parameterization of the initial conditions. 
A more rigorous determination of the uncertainties needed in the lattice EoS is expected within the next few years. 


D. Fluctuations, Freeze-out, and Finite Baryon Density 

As noted, lattice QCD studies have firmly established that at high temperatures and zero baryon chemical potential 
normal hadronic matter turns into a QGP through a smooth, but rapid crossover. However, based on various 
theoretical studies it is generally believed that at large baryon densities and small temperatures the transition from 
hadronic to QGP matter takes place via a first order phase transition. The conjectured point in the temperature- 
baryon chemical potential phase diagram of QCD, at which this first order transition line meets the crossover region is 
known as the QCD critical (end) point [^[^. The QCD critical point is a unique point in the QCD phase diagram 
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beyond which the hadronic and the QGP phase coexist along the first order line. At the QCD critical point, a second 
order phase transition takes place between the hadronic and the QGP phase, resulting in long-range correlations at 
all length scales. 

The large correlation length (^) associated with a nearby critical point manifests itself through increased fluctuations. 
For example, while the second cumulant of conserved charge fluctuations scales as the higher order cubic and quartic 
cumulants grow as and respectively [61]. Furthermore, it has been shown that even qualitative features of 
higher cumulants can signal presence of criticality |62H64] . These cumulants can also be accessed in heavy ion collisions 
via event-by-event fluctuations |65j . In this vein, the search for the QGD critical point in the Relativistic Heavy Ion 
Gollider’s (RHIC) Beam Energy Scan (BES) program is concentrated on measurements of higher order cumulants of 
the fluctuations These higher cumulants of conserved charge fluctuations are also accessible in lattice QCD 

calculations. Among several conserved charges the net electric charge is of special interest as its fluctuations can 
provide a direction comparison between experimental measurements and lattice QCD |69j . 

Although a direct lattice QCD computation at non-zero baryon (/rs), charge (/rg) or strangeness (^s) chemical 
potentials remains difficult due to the infamous sign problem, higher cumulants of fluctuations of these conserved 
charges can be computed on the lattice using the well established method of Taylor expansion [zniiHi]. In this method 
one expands the QCD logarithm of the partition function, In Z, or the pressure, P = —T ln(Z)/P, in a power series of 
the chemical potentials around vanishing values of the chemical potentials. For the electric charge chemical potential 


2^4 


oo 


E 

n=0 



where 


X^iT) 


1 a"lnZ 
VT^ dip^Q/TT 


(17) 


Here, V and T denote the volume and the temperature respectively. The coefficients Xn known as the generalized 
susceptibilities associated with the specific conserved charge. Since these generalized susceptibilities are defined at 
vanishing chemical potentials, lattice QCD simulations can be used to compute them. In order to obtain these 
susceptibilities at pB 0, one can further Taylor expand in a power series of /tb around /is = 0 


(T’, mb) = E ^ ^kn (T) (^) , where Xkn (T) = g 

These generalized susceptibilities are, in turn, related to the fluctuations 
o ^ 1 ,, n . 1 


(18) 


f^B—0 

of the conserved charges 

X?iT,PB) = ^{NQ) , 

X?iT,aB) = ^((SNof) , 


X2 


(19) 


where Nq is the net (positive minus negative) charge and SNq = Nq — {Nq). 

On the other hand, heavy ion experiments measure various cumulants, such as the mean (M), variance (cr), skewness 
(S), kurtosis (k) etc., of the event-by-event distribution of the net charge |66| at a given collision energy {^/s). These 
cumulants are related to the higher order non-Gaussian fluctuations of conserved charge; as an example, for the net 

charge [55] 


Mq (Vs) = {Nq) , 


Sq{V^) 



(Vs) = {{SNq) 


,{SNq) 

(v 


-3 . 


( 20 ) 


Recent experimental advances in measurements of cumulants of charge fluctuations have placed us in a unique 
situation where, for the first time, lattice QCD computations at non-zero temperatures and densities can be directly 
confronted with the results from heavy ion experiments through the use of appropriate volume-independent ratios of 
cumulants of net charge fluctuations m 


Mq {Vs) 
<^Q (Vs) 
Sq (Vs) cr| (Vs) 
Mq{Vs) 


Xi {T.tJ-s) _ ^Q 
x?(r,^B)" 

Xs (Vmb) _ ^Q 
x?(r,^B)" " 


(2Ia) 

(21b) 
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To obtain any information regarding the location of the critical point in the T — hb phase diagram of QCD from 
experimentally measured cumulants of charge fluctuations, it is essential to relate the experimentally tunable param¬ 
eter ^/s to thermodynamic parameters, namely the freeze-out temperature and freeze-out chemical potential Hg. 
Recently, it has been shown m that it is possible to ext ract these thermal freez e-out parameters and fig by 
comparing first principles lattice calculations for [Eq. 21b and [Eq. 21a , directly to their corresponding 


cumulant ratios in heavy ion collisions. The feasibility of such a procedure has been demonstrated in [7M75]. A 
recent example of such a comparison and subsequent determination of the freeze-out parameters are shown in Figure 


11 (see the figure caption for details) [7H] . 




FIG. 11: (Left) The figure shows a comparison between the lattice QCD results |72) for i?® and the STAR data |68) for 
{Sq<j^)/Mq at ^ = 62.4 GeV. The overlap of the experimental results with the lattice QCD calculations provides an upper 
bound on the freeze-out temperature < 155 MeV. (Right) Lattice QCD results [72] for as a function of /is compared 
with the STAR data [68] for Mqjaq in the temperature range = 150(5) MeV. The overlap regions of the experimentally 
measured results with the lattice QCD calculations provide estimates for the freeze-out chemical potential for a given y/s. 
The arrows indicate the corresponding values of fig obtained from the statistical model fits to the experimentally measured 
hadron yields izzl- 


It is evident from Figure [IT| that the large errors on the experimental values {SQaQ)/MQ, allow at present, only an 
upper bound on the freeze-out temperature . Very recently, an alternative procedure for the determination of 
has been demonstrated |78j . This procedure utilizes the fact that the initially colliding nuclei in heavy ion collisions 
are free of net strangeness and conservation of strangeness under strong interaction ensures that the medium created 
during heavy ion collisions is strangeness neutral. By imposing a strangeness neutrality condition for a homogeneous 
thermal medium (ns) (ubiRs) = 0, the strangeness chemical potential, /is, can be determined by performing a Taylor 
expansion of the net strangeness density, (ns) (/is,/is) [ZH] 


Mb 


si(T) + s3(r)(^)' + o 



( 22 ) 


Since the coefficients Si, S 3 , etc. consist of various generalized baryon, charge and strangeness susceptibilities defined 
at vanishing chemical potentials, they can also be obtained from standard lattice QCD computations at zero chemical 
potentials. The leading order si(T) coefficient for /is//iB is shown in Figure [l^ (left). It is interesting that com¬ 
parisons of these lattice results with the predictions from the hadron resonance gas model reveal that the inclusion 
of only experimentally observed hadrons fails to reproduce the lattice data around the crossover region. However, 
the inclusion of additional, unobserved strange hadrons predicted within the quark model provides a much better 
agreement with lattice results, hinting that these additional hadrons become thermodynamically relevant close to the 
crossover temperature m- Other lattice thermodynamics studies also indicate that additional, unobserved charm 
hadrons also become thermodynamically relevant close to the QCD crossover m- 
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In heavy ion collisions, the measured relative yields of the strangeness S antibaryons to baryons at the freeze-out 
are determined by the thermal freeze-out parameters {T-^, [77] 


Rh{Vs) = exp 




(23) 


Thus, by fitting the experimentally measured values of Ra, Rb. and Rq, corresponding to [S'! = 1,2 and 3, at a given 
it is easy to determine the matching these experimentally extracted values of /i 5 //i 5 with 

lattice QCD results for fig/^J^B as a function of temperature, one can determine the freeze-out temperature . Figure 
12 (right) demonstrates this procedure. Not surprisingly, the inclusion of additional unobserved strange hadrons in the 
hadron resonance gas model leads to very similar values of the freeze-out temperatures as obtained using the lattice 
data. However, including only the hadrons listed in the Particle Data Group tables |HT| yields freeze-out temperatures 
that are 5 — 8 MeV smaller. 
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FIG. 12: (Left) Lattice QCD results for /rs/ps at the leading order, i.e. Si(T) (see Eq. (^). The dotted line (PDG-HRG) 
shows the results of hadron resonance gas model containing only hadrons listed by the Particle Data Group. The solid line 
(QM-HRG) depicts the result for a hadron gas when additional, yet unobserved, quark model predicted strange hadrons are 
included [78]. The shaded region indicate the chiral crossover region Tpc = 154(9) MeV. (Right) The figure shows a comparison 
between the experimentally extracted values of {/ig/fig, Hg/T^) (filled points) with the lattice QCD results for fis/^I'B (shaded 
bands) [78]. The lattice QCD results are shown for fis/T = fig/TR The temperature range where lattice QCD results match 
with fig/fig provide the values of T-^, i.e. T = 155(5) MeV and 145(2) MeV for ^/s = 39 GeV and 17.3 GeV, respectively. 


Lattice QCD calculations can also be used to locate and establish the existence of the QCD critical point. Early 
results are based on calculations on rather small and coarse lattice using only a l-link standard staggered fermion 
discretization scheme [80] . It has been pointed out that the method used in this calculation, the determination of 
Lee-Yang zeroes, also suffers from an overlap problem rather than a sign problem and may lead to spurious signatures 
for a critical point |82j . Calculations using a formulation of hnite-density QCD with an imaginary chemical potential, 
also performed on lattices with only four sites in the temporal direction, do not hnd any evidence for the existence 
of a critical point |H3]. The most systematic searches for a critical point at present are based on the Taylor series 
expansion of the QCD partition function |84l|85|. For vanishing electric charge and strangeness chemical potential 
one can expand the pressure P in terms of {fiB/TY', the expansion coefficients being cumulants of net-baryon number 
fluctuations, i.e. generalized baryon number susceptibilities Xni"R) 


P{T,flB) 


OO . 

n—0 


{T){ 


MB V 
T / 


(24) 


This series expansion has a radius of convergence which may be estimated using a finite, typically small, set of values 
for the generalized baryon number susceptibilities. Subsequent estimators 


r-n 



/ n(n- I)xg 
Xu+2 


n even , 


(25) 


may stay finite or diverge in the limit n —>■ oo. Current estimates for the radius of convergence |86j . also based on 
calculations with the unimproved I-link staggered fermion action and moderately light quark masses (M.^. ~ 230 MeV), 
suggest for the coordinates of the critical point {T^/Tp,., fig/T^) = (0.94(1), 1.8(1)). However, true systematic errors 
for these calculations are difficult to estimate, and this topic is currently under active investigation. 
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IV. TRANSPORT PROPERTIES AND HEAVY QUARKS 

This review focuses on recent results on bulk QCD thermodynamics with physical quark masses and on fluctuations 
of conserved charges. These two topics are of immediate importance for the ongoing experimental studies on the 
phase structure of strong-interaction matter. However, a comprehensive understanding of the strongly coupled nature 
of QGP for temperatures Tpc ^ T < 2Tpc requires lattice QCD calculations of its color screening and transport 
properties. Many of these calculations are performed in the quenched approximation, in which the dynamical fermion 
loops are neglected. Furthermore, full continuum extrapolations have not yet been performed, but they are often 
performed on large lattices that are close enough to the continuum to yield meaningful results. We give a brief 
summary of these important calculations and the insights that they provide. 


A. Color Screening 

Matsui and Satz pointed out that the force between heavy quarks inside a QGP is screened due to presence of color 
charges and eventually leads to the dissolution of quarkonia, bound states of heavy quark and anti-quark, such as J/i/j, 
rjc, T etc. [57]. Lattice QCD calculations of the potential between two infinitely heavy static quarks have established 
this color screening mechanism |55J |5S|. However, quantitative understanding of the dissociation temperatures of 
quarkonia require knowledge of the spectral functions of quarkonia. Extraction of these real (Minkowski) time quanti¬ 
ties from the Euclidean time quarkonia correlation functions measured on the lattice require analytic continuations to 
Minkowski time. Since the number of lattice points along the Euclidean time direction are limited, such analytic con¬ 
tinuations are usually performed using a Bayesian method, such as the Maximum Entropy Method uniiii]- Reliable 
analytic continuation via the Maximum Entropy Method demands lattice data at large numbers of Euclidean time 
points, and, hence, lattices with large temporal extents. Galculations of charmonium spectral functions in quenched 
QGD have been performed on large lattices, close to the continuum limit [55]. These calculations suggest that both P- 
and S-wave ground state charmonia disappear in a QGP at temperatures T > l.STpc- Calculations with 2-1-1 flavors 
of dynamical quarks, but with un-physically heavy pion masses, have started |93L I94j and lead to similar conclusions. 
Lattice QCD study of the spatial correlation functions of charmonia with 2-1-1 dynamical flavors having almost 
physical quark masses has also provided indirect evidence that the ground state charmonia cease to exist inside QGP 
for T > l.STpc 1^. At present, lattice spacings (a) used in the state-of-the-art lattice QCD calculations at non-zero 
temperatures are still too large such that the bottom quark mass (m/,) in lattice units are ami, ^ 1- This leads to 
large cut-off effects for studies related to bottomonia spectral functions. Thus, currently one uses a hybrid approach 
where the heavy bottom quarks are treated within the non-relativistic QCD (NRQCD) approximation. However, the 
NRQCD formalism does not possess a proper continuum limit. These calculations inaiMiinT] suggest that the ground 
state S-wave bottomonium survive up to ^ 2Tpc, but the implications for the P-wave states are still unclear. 

B. Transport Coefficients 

Calculations of transport coefficients of QGP using Euclidean time lattice QCD also require analytic continuation to 
real Minkowski time, and therefore also depend on the extraction of spectral functions. Shear and bulk viscosity can 
be obtained from Euclidean time correlation functions of the energy-momentum tensor. Since the energy-momentum 
tensor operator is dominated by purely gluonic correlation functions these correlation functions tend to be very noisy 
and calculations of shear and bulk viscosities on the lattice remain extremely challenging. Although, there were 
attempts to calculate shear and bulk viscosities for a pure gauge theory [55lH01j , so far, there is no lattice calculation 
of these quantities for realistic QCD. 

Transport coefficients determined solely by quark operators, such as the electrical conductivity, are more accessible 
from present day lattice QCD. Calculations of electrical conductivity (cr) within the quenched approximation are quite 
advanced and results close to the continuum limit exist m- A recent calculation |103j at three temperatures in the 
range T = f.lT'pc — 1.5Tpc constrains the value to a rather narrow range, 0.2Cem ^ cr/T < 0.4Cem, where Cem denotes 
the sum over the squared electric charges of quarks. More realistic calculations of tr with dynamical fermions have also 
started to become available in recent years |104H10f)| and show a striking drop in the value of cr/T when a approaching 
Tpc from high temperatures [10511106] . Similar calculations also provide the charge diffusion constant [106] and the 
thermal di-lepton production rate in QGP [10211103] . 

Spectral functions associated with the vector current for the charm quark also provide access to the charm quark 
diffusion constant (D) [S^]- The momentum diffusion coefficient of an infinitely heavy quark, k = 2T^/D, can 
also be extracted from the correlation function of purely gluonic operators under the heavy quark effective theory 
approximation m- Current results [10811109] on the momentum diffusion constant of infinitely heavy quark are 
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consistent, yielding k ~ 2.5T^. The corresponding diffusion constant D ~ 0.8/T is about a factor two larger than the 
charm quark diffusion constant extracted using the charm vector current correlation function |92j . 

The jet quenching parameter q, an important ingredient in the analysis of energy loss of jets traversing QGP, may 
also become accessible to lattice QCD calculations. Calculation of this quantity involves correlation functions of Wilson 
lines along the light-cone and extraction of that from Euclidean lattice QCD again demands analytic continuation to 
real time. One way of performing such analytic continuation can be justified at very high temperature where the weak 
coupling expansion is valid |110j . Following this proposal, a part of the non-perturbative contributions to the q was 
calculated recently within the dimensionally reduced effective theory, electrostatic QCD |111) . If this proposal can be 
extended to full QCD and down to the truly non-perturbative regime close to the QCD crossover then it will open up 
a new avenue for lattice QCD calculations that will directly impact the phenomenology of strongly interacting matter 
probed in heavy ion collisions. 


V. CONCLUSIONS AND OUTLOOK 

Calculations of the fundamental thermodynamic quantities of QCD as presented in this review have reached a 
significant milestone. The crossover temperature and equation of state have been calculated with physical values 
for the light and strange quark masses and separate, reliable continuum extrapolations have been performed by two 
collaborations. The transition is firmly established as a crossover, as evident in both the equation of state results and 
observation that chiral susceptibility is independent of volume. The latter result has been achieved with fermions 
that differ in their approach to chiral symmetry: a staggered fermion action in which full chiral symmetry is restored 
in the continuum, and the domain wall action in which chiral symmetry is preserved for finite lattice spacings. 
Because the transition is a crossover, the definition of a transition temperature is quantity dependent. However, 
the chiral condensate, which is the order parameter for the phase transition in the chiral limit, is a natural choice. 
For this quantity there is also remarkable agreement among the recent calculations despite significant differences in 
the analysis methods used. A fit to the inflection point in the renormalized chiral condensate with the stout action 
yielded Tp^ = 155 ± 2 ± 3 MeV. A more sophisticated analysis involving scaling fits to the 0{N) universality class for 
a constrained continuum extrapolation to the HISQ and asqtad actions produced a value of Tpc = 154 ± 8 ± 1 MeV. 
The general range for Tpc has also been reproduced in a domain wall calculation with physical quark masses and a 
lattice spacing Nr = 8. 

The close agreement between different actions and analysis methods extends to the equation of state, where con¬ 
tinuum extrapolations with the stout and HISQ actions agree to within their respective errors over the temperature 
range 130-400 MeV. Although a small difference begins to develop at the higher end of this temperature range, it is 
not yet known whether this will lead to a more significant difference above this temperature range. At this time the 
overall precision and estimated accuracy of the crossover temperature and equation of state at zero baryon density 
appear to be sufficient to meet the needs of the heavy ion community. The need for future improvements from the 
lattice will depend upon the sophistication of the phenomenological tools and the desired accuracy for extracted 
physics parameters. 

Calculations of fluctuations on the lattice are a more recent development, and this area has received considerable 
attention and resources only within the past few years. The ability to calculate freeze-out curves and net-charge 
moments that can be directly compared with heavy ion experiments represents a significant advance, and one the 
will hopefully elucidate the location and signatures of the critical point. At this time predictions of the critical point 
in the T-fj, plane are highly uncertain, and significant advances are required in both computing power and algorithm 
efficiency if this goal is to be attained within the next several years. 

Finally, calculations of light and heavy quark bound states, diffusion, and now jet quenching hold considerable 
promise for the future, and one can expect significant results to follow when full QCD calculations with physical 
quark masses are possible. 

This work is supported in part through contracts No. DE-AC52-07NA27344 and No. DE-SC0012704 with the U.S. 
Department of Energy and NSF Crant No. PHYlO-034278. 
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